Bulletin of Mathematical Biology

Bladder cancer presents significant clinical challenges due to its complex immune microenvironment and highly heterogeneous response to treatments. To create accurate, individualized models of disease progression, we first construct a system of Ordinary Differential Equations (ODEs) that captures tumor-immune interactions. We address the challenge of estimating unknown parameters by performing a rigorous comparative analysis of two heuristic optimization methods: Differential Evolution (DE), a robust global optimization algorithm, and Physics-Informed Neural Networks (PINN), a novel machine learning framework that embeds ODE constraints into its loss function. Our findings provide a critical evaluation of the computational efficiency and accuracy of each method for parameterizing biological ODE systems. This study validates the power of hybrid machine learning approaches in mathematical oncology, yielding a data-driven model of bladder cancer progression with direct potential for optimizing personalized treatment strategies. Author summaryBladder cancer remains a major global health threat, characterized by highly unpredictable responses to treatment and a high likelihood of recurrence. To better predict how a patients disease will progress, researchers use mathematical models that simulate the "war" between cancer cells and the immune system. However, these models are only useful if they can be accurately tuned to a specific patients data--a process called parameter estimation. This task is notoriously difficult because clinical data is often sparse and noisy, making it hard to find the right settings for the model. In this study, we developed a novel computational framework that combines a traditional "survival of the fittest" optimization algorithm (Differential Evolution) with Physics-Informed Neural Networks (PINNs), a specialized architecture designed to embed physical constraints directly into the learning process. By "teaching" the AI the underlying biological laws of cancer growth, our hybrid approach can accurately estimate a patients unique disease parameters even when raw data is limited. We validated this method using a "virtual patient" system derived from real-world clinical trials. Our results show that this hybrid approach provides a more robust and reliable way to personalize cancer models, offering a powerful new tool for doctors to simulate and optimize individual treatment plans before they are even administered.

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

Cancer-associated fibroblasts (CAFs) are a central component of the tumor microenvironment that facilitate a supportive niche for cancer progression and metastasis. Experimental evidence suggests that CAFs can facilitate estrogen-independent tumor growth, thereby reducing the efficacy of anti-hormonal therapies. Understanding and quantifying the complex interactions between tumor cells, hormonal signalling, and the microenvironment are crucial for designing more effective and individualized treatment strategies. We propose a mathematical framework to explore the influence of CAFs on ER+ breast cancer progression and to evaluate strategies to mitigate their impact. We develop a system of nonlinear ordinary differential equations that substantiates the experimental observations by providing a mechanistic basis for the role of CAFs in regulating estrogen-independent growth dynamics. We then employ optimal control theory to evaluate distinct therapeutic approaches involving monotherapy or combinations of: (i) inhibition of tumor-to-CAF signaling, (ii) inhibition of CAF-to-tumor proliferative signaling, and (iii) endocrine therapy. Taken together, our results demonstrate that CAF-targeted strategies can enhance treatment efficacy across various estrogen dosing regimens. Our study provides new insights into the potential of CAF as a therapeutic target that could help to improve existing approaches for endocrine therapies.

Biological tissues grow at rates that depend on the geometry of the supporting tissue substrate. In this study, we present a novel discrete mathematical model for simulating biological tissue growth in a range of geometries. The discrete model is deterministic and tracks the evolution of the tissue interface by representing it as a chain of individual cells that interact mechanically and simultaneously generate new tissue material. To describe the collective behaviour of cells, we derive a continuum limit description of the discrete model leading to a reaction-diffusion partial differential equation governing the evolution of cell density along the evolving interface. In the continuum limit, the mechanical properties of discrete cells are directly linked to their collective diffusivity, and spatial constraints introduce curvature dependence that is not explicitly incorporated in the discrete model. Numerical simulations of both the discrete and continuum models reproduce the smoothing behaviour observed experimentally with minimal discrepancies between the models. The discrete model offers further individual-level details, including cell trajectory data, for any restoring force law and initial geometry. Where applicable, we discuss how the discrete model and its continuum description can be used to interpret existing experimental observations.

Vaccine-acquired immunity plays an important role in controlling the spread of many infectious diseases; however, vaccine efficacy can diminish over time. This work uses a mathematical model to study the effects of waning vaccination-acquired immunity on infection incidence. With an SEIR-type compartmental model that considers both vaccinated and unvaccinated populations (and their mixing), we present mathematical conditions under which vaccinated individuals drive ongoing growth in infections, i.e., over half of the daily incidence arises from vaccinated individuals. Analysis of a mathematical model of COVID-19 spread in the state of Colorado suggests how and for what duration vaccinated individuals could have sustained such growth. Importantly, our model demonstrates that, despite potential for brief vaccinated-driven periods of growth in infections, which occur among unvaccinated-driven periods of growth in infections, increased vaccination coverage always reduces total cases and total hospitalizations. This work provides insight into how waning immunity in vaccinated populations can contribute to ongoing infection incidence and demonstrates the value of complementary interventions to prevent disease spread in vaccinated populations.

We present an approach to analysing cell homeostasis using a bond graph modelling approach that ensures that the conservation laws of physics (conservation of mass, charge, and energy, respectively) are satisfied for the interdependent biochemical, electrical, mechanical, and thermal energy storage mechanisms operating within the cell. We apply the bond graph approach to several cell membrane transport mechanisms and then consider how physics constrains intracellular electrolyte homeostasis for enterocytes (the epithelial absorptive cells of the gut). The model includes the electrogenic sodium-potassium ATPase pump (NKA), the glucose transporter (GLUT2), and an inwardly rectifying potassium channel, all in the basolateral membrane, and the electrogenic sodium-driven glucose transporter (SGLT1) in the apical membrane. Glycolysis converts the imported glucose to ATP to drive NKA. For specified levels of sodium, potassium, and glucose in the blood, the model demonstrates how enterocytes absorb sodium and glucose from the gut and transfer glucose to the blood while maintaining the membrane potential and homeostasis of intracellular sodium and potassium. The Gibbs free energy available from the ATP hydrolysis ensures that the cell operates as a sodium battery with a high external to internal ratio of sodium concentration in order to provide the energy for many other cellular transport processes. We show that the 3:2 stoichiometry of Na+/K+ exchange in NKA, coupled with 2:1 Na+/glucose cotransport in SGLT1, a 1:2:2 ratio between glucose consumption and ATP and water production in glycolysis, and K+ and glucose efflux through Kir and GLUT2, respectively, provides a balanced system that maintains homeostasis of intracellular Na+, K+, glucose, ATP and water, and homeostasis of the membrane potential, under varying levels of transport of glucose from the gut to the blood. We also show how the flux expressions for SLC transporters, ATPase pumps and ion channels can all be expressed in a consistent and thermodynamically valid way.

We introduce a mechanistic, nonlocal tumour-growth model designed specifically to capture explosive dynamics that are not adequately explained by standard logistic reaction-diffusion descriptions. The motivation is empirical: the universal scaling law reported in [1] provides compelling cross-sectional evidence of superlinear tumour activity versus tumour burden, but as a phenomenological relationship it does not by itself supply a dynamical mechanism, nor does it rigorously describe how explosive growth emerges, how fast it develops, or how spatial interactions and tissue boundaries influence it. Our model addresses this gap by incorporating nonlocal proliferative feedback--cells respond to a spatially aggregated neighbourhood signal--and a singular, Kawarada-type acceleration that produces "quenching": tumour density stays bounded while the proliferative drive becomes unbounded as the aggregated signal approaches a critical threshold. This offers a concrete mechanistic route to explosive escalation consistent with physical boundedness. We analyse the model under no-flux (Neumann) boundary conditions, appropriate for reflecting tissue interfaces. In the spatially homogeneous setting we prove finite-time onset of the explosive regime and obtain explicit rates for how rapidly it is approached. For spatially heterogeneous perturbations we derive a transparent spectral stability theory showing how the interaction kernel selects spatial scales and how the singular acceleration tightens stability margins as the explosive threshold is approached. These results provide interpretable links between nonlocal interaction structure, boundary effects, and the emergence of rapid growth. Finally, to connect mechanism to data in the spirit of [1], we embed the model in a Bayesian inference framework that treats the interaction kernel and the acceleration strength as unknown and learned from tumour-growth observations. This enables uncertainty-aware estimation of explosive onset times, escalation rates, and stability margins, while positioning the scaling law of [1] as an observable signature that our mechanistic model can explain and quantify rather than merely fit.

Variability in gene expression among single cells and growing cell populations can arise from the stochastic nature of protein synthesis, which often occurs in random bursts. This study investigates the variability in the expression of a growth-sustaining protein, whose concentration is regulated by a negative feedback loop due to cell growth-induced dilution. We model the distribution of protein concentration using a Chapman-Kolmogorov equation for single cells and a population balance equation for growing cell populations. For single cells, we derive an explicit solution for the protein concentration distribution in state space and represent it as a Bessel function in Laplace space. For growing populations, we find that the distribution satisfies a Heun differential equation with singular boundary conditions. By addressing the central connection problem for the Heun equation, we quantify the population-level protein distribution and determine the Mathusian parameter, which characterizes population growth. This work provides a comprehensive analytical framework for understanding how stochastic protein synthesis impacts gene expression variability and population dynamics.

In this research, we create a new fractional-order SEIHRD framework to examine how the Nipah virus moves from one species to another (zoonotic spillover) and how it later spreads throughout a community (via contact with one another) or in a hospital or isolation situation (via entering into a hospital or being placed under quarantine). We used the fractional-derivative formulation of the SEIHRD model to demonstrate memory-based effects related to the progression of an infection and also reflect time-distributed effects associated with surveillance and control measures placed on an infected patient. We first demonstrated that the basic epidemiologic properties of the model were consistent by showing that the solutions of the SEIHRD differential equations will always yield positive and bounded solutions within biologically relevant parameter ranges. We then established the well-posedness of this model by transforming the SEIHRD differential equations into an equivalent integral operator and applying various fixed-point arguments to demonstrate that there will always be unique solution(s) to the SEIHRD differential equations. To evaluate the threshold parameter for the transmission of Nipah virus within a given population we calculated the threshold level through the next generation method to determine the expected number of secondary infections from a new or chronically infected host. One of the main contributions of this work is to include an analysis of the robustness of a given solution to all potential perturbations (i.e., Ulam-Hyers and generalized Ulam-Hyers stability). In addition, we provide analytic results guaranteeing that small perturbations due to approximate modeling, numerical approximation (discretization), or the lack of data fidelity will produce controlled deviations in the solutions. To finish this project, we perform a global sensitivity analysis on uncertain coefficients to evaluate their contribution to the uncertainty of each coefficient and to find out the coefficients that most strongly influence major outcome metrics. This will allow us to develop a priority order for prioritizing spillover control (reduction of human contact and/or isolation), contact reduction, and expenditure of resources towards isolation-related interventions. The resulting framework converts fractional epidemic modeling from a descriptive simulation to a replicable method with robustly defined behavior and equal response prediction.

Cytokine-mediated communication is a central mechanism by which immune cells coordinate activation, differentiation and proliferation. While mechanistic reaction-diffusion models provide detailed descriptions of cytokine secretion and uptake at the cellular scale, their computational cost limits their applicability to large and densely packed cell populations. Previously employed approximations of cytokine diffusion fields rely on assumptions that neglect the influence of cellular geometry and volume exclusion. In this work, we study a macroscopic description of cytokine diffusion and reaction dynamics based on homogenization techniques, rigorously linking microscopic reaction-diffusion formulations to effective continuum models. The resulting homogenized equations replace discrete responder cells with a continuous density, while retaining essential features of cellular uptake and excluded-volume effects. Further, we show that in regimes with approximate radial symmetry, classical Yukawa-type solutions emerge as limiting cases of the homogenized model, provided appropriate correction factors are included. Overall, our approach allows efficient multiscale modeling of cytokine signaling in complex immune-cell environments.

Neuronal activity alters extracellular ion concentrations and electric potentials. Ephaptic effects refer to the feedback influence that these extracellular changes can have on neuronal activity. While electric ephaptic effects occur on a fast timescale due to extracellular potential perturbations, ionic ephaptic effects are driven by slower, accumulative changes in ion concentrations. Among the previous computational studies of ephaptic effects, the vast majority have focused exclusively on electric effects, while ionic ephaptic effects have largely been neglected. In this work, we present an electrodiffusive computational framework consisting of two-compartment neurons that interact via a shared extracellular space. By accounting for both electric potentials and ion-concentration dynamics in a self-consistent manner, our framework enables us to explore the relative roles of electric and ionic ephaptic effects. Through numerical experiments, we demonstrate that ionic and electric ephaptic interactions play very different roles. While ionic ephaptic interactions increase population firing rates, electric ephaptic interactions primarily drive subtle shifts in spike timing. Furthermore, we show that these spike shifts cause the phase difference (the distance in spike times between a small collection of neurons) to converge to a stable, unique phase difference, which we coin the ephaptic intrinsic phase preference. Author summaryNeurons predominantly communicate through synapses: specialized contact points where a brief electrical signal, known as a spike or action potential, in one neuron influences another. Neurons generate these spikes by exchanging ions with the surrounding extracellular space. This way, spiking neurons alter extracellular ion concentrations and electric potentials. Since neurons are sensitive to such changes in their environment, they can also influence one another indirectly through the shared extracellular medium. This form of non-synaptic interaction is known as ephaptic coupling. Most computational models of neuronal activity neglect ephaptic interactions, and those that include them typically consider only electric effects while ignoring ionic contributions. As a result, the relative roles of electric and ionic ephaptic effects remain poorly understood. Here, we introduce a computational framework that accounts for both mechanisms in a self-consistent way. Our results show a functional distinction: ionic ephaptic effects act slowly, regulating population firing rates, whereas electric ephaptic effects act on millisecond timescales and subtly shift spike timing. These shifts cause spike-time differences between neurons to converge to a stable value, a phenomenon we call ephaptic intrinsic phase preference.

Evolutionary cancer therapy (ECT) applies principles of evolutionary game theory to prolong the effectiveness of cancer treatment by curbing the development of treatment resistance. It was shown to increase time to progression while decreasing the cumulative drug dose. ECT individually tailors treatment schedules for patients based on their cancer dynamics and, thus, requires regular follow-up and precise measurements of the cancer burden. The current literature on ECT often overlooks clinical realities, such as rather long intervals between tests, possible appointment delays and measurement errors, in the development of the treatment protocols. In this study, we assess the clinical feasibility of ECT for metastatic non-small cell lung cancer (NSCLC). We create virtual patients with cancer dynamics described by the polymorphic Gompertzian model, based on data from the START-TKI clinical trial. We assess the effects of longer test intervals, measurement error and appointment delays on the expected time to progression under the evolutionary therapy protocols. We show that a higher containment level, although it increases time to progression in the models predictions, may lead to premature treatment failure in the presence of measurement error and appointment delay. Further, we show that the ECT protocol with a single containment bound is more robust to the clinical realities than the protocol with two bounds. Finally, we show that a dynamically adjusted treatment protocol can be beneficial for individual patients, but requires a thorough follow-up. This study contributes to the design of a clinical trial and the future clinical implementation of evolutionary therapy for NSCLC.

We present the Toroidal Search Algorithm (TSA), a novel population-based metaheuristic optimization method inspired by the topology of a torus. Conventional metaheuristics frequently suffer from boundary stagnation, a phenomenon that severely degrades performance in bounded and high-dimensional search spaces. TSA addresses this limitation by embedding the search domain into a toroidal geometry, thereby eliminating artificial boundaries and enabling continuous cyclic exploration. Beyond boundary handling, TSA uses winding numbers to capture the history of agent movement across periodic dimensions, which are exploited to adaptively refine local search. A modified sigmoid control function regulates the transition between global and local search. Performance of TSA is evaluated on a collection of unimodal and multimodal benchmark functions at various dimensions. It consistently outperforms established metaheuristics. Notably, TSA demonstrates exceptional robustness to increasing dimensionality, maintaining fast convergence and low variance where competing methods deteriorate. To assess real-world applicability, we apply TSA to an inverse problem from mathematical oncology. With both synthetic and clinical data, TSA reliably recovers physiologically plausible parameters with greater stability and predictive accuracy than competing algorithms. These results demonstrate that TSA is a powerful and robust tool for large-scale global optimization in computational modelling applications. Striking Image O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=200 SRC="FIGDIR/small/709766v1_ufig1.gif" ALT="Figure 1"> View larger version (131K): org.highwire.dtl.DTLVardef@1b2c994org.highwire.dtl.DTLVardef@d045a8org.highwire.dtl.DTLVardef@18d296corg.highwire.dtl.DTLVardef@9a972d_HPS_FORMAT_FIGEXP M_FIG C_FIG Image generated with Google Gemini.

Cerebrospinal fluid (CSF) is a Newtonian fluid that bathes the brain and spinal cord and oscillates in response to the physiological periodic changes in brain volume, of which the cardiac cycle is a major driver. Understanding this motion is essential for clarifying its contribution to solute transport, waste clearance, and drug delivery. In this work, we study oscillatory and steady streaming flow in the cranial subarachnoid space using a lubrication-based theoretical framework. The model represents the cranial CSF compartment as a thin fluid layer bounded internally by the brain surface and externally by the dura, driven by time-dependent brain surface displacements. We first derive simplified governing equations for flow over an arbitrary smooth sphere-like brain surface and obtain analytical solutions for an idealised spherical geometry with uniform displacements. We then incorporate realistic displacement fields reconstructed from MRI measurements in healthy subjects and solve the reduced equations numerically. The results show that oscillatory forcing produces a steady streaming component that may enhance solute transport compared with diffusion alone. This work provides a mechanistic description of the flow generated by physiological brain motion and highlights the potential presence of steady streaming in cranial subarachnoid fluid dynamics.

Mathematically integrating genetics, development, and evolution is a longstanding challenge. Here I develop general mathematical theory that integrates sexual, discrete, multilocus genetics, development, and evolution. This yields an exact method to describe the evolutionary dynamics of allele frequencies and linkage disequilibria in multilocus systems and the associated evolutionary dynamics of mean phenotypes constructed via arbitrarily complex developmental processes. The theory shows how development affects evolution under realistic genetics, namely by shaping the fitness landscape of allele frequencies and linkage disequilibria and by constraining adaptation to an admissible evolutionary manifold (high dimensional region on the landscape) where mean phenotypes, phenotype (co-)variances, and higher moments can be developed. I derive a first-order approximation of this exact method, which yields equations in gradient form describing change in allele frequency, linkage disequilibria, and mean phenotypes as constrained, sometimes-adaptive topographies. Both the exact and approximated equations describe long-term phenotypic and genetic evolution, including the evolution of mean phenotypes, phenotype covariance matrices, "mechanistic" additive genetic cross-covariance matrices, and higher moments. I provide worked examples to illustrate the methods. The theory obtained is referred to as evo-devo dynamics, which can be interpreted as an extension of population genetics, with some similarities to quantitative genetics but with fundamental differences. The theory provides tools to re-assess empirical observations that have been paradoxical under previous theory, such as the maintenance of genetic variation, the paradox of stasis, the paradox of predictability, and the rarity of stabilising selection, which appear less paradoxical in this theory.

Haptotaxis is an understudied form of directed cell migration in which movements are biased by gradients of immobilized ligands. For example, fibroblasts and other mesenchymal cells sense and respond to gradients of extracellular matrix (ECM) composition, which is relevant during tissue morphogenesis and repair. As a step towards understanding how haptotactic gradients spatially bias cell adhesion, intracellular signal transduction, and cytoskeletal dynamics, we formulated a phase field model of whole-cell migration, in which the occupancy of potential adhesion sites changes stochastically with time. With careful assignment of parameter values, the model predicts significant haptotactic bias for adhesion-site gradient steepness of a few percent across the cell. We then used the model to predict how the cells removal of surface-bound ECM ligand (as observed in experiment) and/or the presence of a competing, chemotactic gradient influence(s) haptotactic fidelity. An emergent principle is that gains in directional persistence naturally offset losses of directional bias, at the cost of greater cell-to-cell heterogeneity of the response. In the case of orthogonally oriented gradients, this offset manifests as a remarkable robustness of the multi-cue response.

Cerebrospinal fluid (CSF) circulates around and through the brain, supporting neural homeostasis by regulating the extracellular chemical environment. Yet the physical mechanisms governing CSF-driven solute transport remain poorly understood, limiting the design of diagnostic and therapeutic strategies targeting brain clearance and drug delivery. Pulsatile CSF flow in the cranial subarachnoid space (cSAS), is driven by cardiac, respiratory, and sleep-related vasomotion. Over longer timescales weaker steady flows, such as inertial steady streaming, Stokes drift, and production-drainage flow, may contribute to solute transport, but their role and relative importance remain unclear. Here, we develop a simplified two-dimensional model of CSF flow and solute transport in the cSAS using lubrication theory. Through multiple-timescale and asymptotic analyses, we derive a reduced long-time transport equation in which advection is governed by the Lagrangian mean velocity, incorporating steady streaming, production-drainage flow, and Stokes drift. Analysing three physiologically relevant case studies, we show that steady flows can substantially reshape concentration profiles, enhance dispersion, and alter clearance efficiency. Our results clarify the mechanisms underlying CSF-mediated transport, predict distinct regimes in humans and mice, and highlight the importance of subject-specific physiological parameters when interpreting contrast-agent and intrathecal drug-delivery studies.

The Bayesian Ordered Lattice Design (BOLD) method for Phase I clinical trials is extended to address an important challenge. It is widely understood that conventional Phase I trial designs are not consistently effective in determining safe and active dose levels. The US FDA launched the Project Optimus, aimed at reforming the paradigms of dose optimization and selection. We propose a backfill BOLD design (BF-BOLD) that centers on BOLD for dose-finding but also adds an activity evaluation for each patient. Our method for determining the optimal biological dose (OBD) first involves identifying the maximum tolerated dose (MTD) and then assessing activity rates among dose levels below the identified MTD. This approach is straightforward and does not require complex statistical modeling. The results of the simulation indicate that performing dose-finding trials with backfilling can both enhance safety and activity assessment, thereby improving treatment sustainability while also preserving the potential for efficacy of the Recommended Phase II Dose (RP2D). We also demonstrate the applicability of the backfill design for reducing overdose rates, and as a more attractive alternative to small-scale dose expansion trials that follow dose escalation. Backfill designs are an important design approach for early phase trials.

Mass-conserved reaction-diffusion systems are used as mathematical models for various phenomena such as cell polarity. Numerical simulations of this system present transient dynamics in which multiple stripe patterns converge to spatially monotonic patterns. Previous studies indicated that the transient dynamics are driven by a mass conservation law and by variations in the amount of substance contained in each pattern, which we refer to as "pattern flux". However, it is challenging to mathematically investigate these pattern dynamics. In this study, we introduce a reaction-diffusion compartment model to investigate the pattern dynamics in view of the conservation law and the pattern flux. This model is defined on multiple intervals (compartments), and diffusive couplings are imposed on each boundary of the compartments. Corresponding to the transient dynamics in the original system, we consider the dynamics around stripe patterns in the compartment model. We derive ordinary differential equations describing the pattern dynamics of the compartment model and analyze the existence and stability of equilibria for the reduced ODE with respect to the boundary parameters. For a specific parameter setting, we obtained results consistent with previous studies. Moreover, we present that the stripe patterns in the compartment model are potentially stabilized by changing the parameter, which is not observed in the original system. We expect that the methodology developed in this paper is extendable to various directions, such as membrane-induced pattern control.

We extend previous models of cumulative cultural evolution by incorporating structured populations and social networks. We examine how connectivity and network topology shape the accumulation of cultural complexity under unbiased (copy randomly), indirectly biased (copy successful individuals), and directly biased (copy successful traits) transmission. We consider random, scalefree, and small-world networks, as well as the communication structures introduced by Mason and Watts, and derive analytical approximations for the homogeneous case. We find that the effects of social structure depend strongly on the form of transmission bias. Under unbiased transmission, network effects are weak except at very low connectivity. Under indirect bias, cultural complexity increases with connectivity, whereas direct bias shows optimal performance at intermediate connectivity, reflecting a trade-off between diffusion and diversity. Differences across topologies are generally modest once the average degree is fixed. Overall, our results show that no single social structure universally promotes cumulative cultural evolution; instead, its effects depend primarily on the dynamics of learning and innovation.